File a8_ dissimilarities and clustering By Dr. Carr Purpose Introductions Dissimilarities and distance options Hierarchical clustering Dendrograms Sections 1. Computing Dissimilarities 1.1 Euclidean distance 1.2 Similarities, transformations and correlation 1.3 1 - cos(angle) between two vectors 1.4 Mutual information 1.5 Manhattan (city block) distance points on a grid 1.6 Jaccard coefficient for binary data 1.7 BLended APples and ORanges dissimilarity 1.8 Assessing the dissimilarity of different 1.9 Dissimilarities and nonlinear dimension reduction 2. The R distance function dist() 3. Four general approaches to clustering 4. Stepping thru an agglomeration clustering example and rank ordering of multivariate data 4.1 The hclust() options and output 4.2 Stepping through single link clustering 4.3 Merge order for dendrograms and linear orderings of multivariate cases 5. Dendrograms and cluster tree variations 6. A dendrogram for a subset of gene expression data 7. Defining clusters for the gene expression data 8. Heat Maps or Color Matrices - Cluster visualization link 9. Clustering different data sets 9.1 iris data plot 9.2 heptathlon pairs plot 9.3 Hypersurface data 9.3.1 Hypersurface generation 9.3.2 clustering and a scatterplot matrix 9.3.3 Rotating ray plot view with cluster means 10. K means Due: 5. One plclust example 6.6 The hypertree dendrogram 8.1 Heat map example 9.1 The overplotted plclust plot of the iris data 9.2 The pairs plot of heptathon data 9.3.2 The pairs plot of the hypersurface 9.3.3 A single static plot of the ray glyphs 10. The pairs plot Obtain From the class web site hypertree.r gemIdGene.csv # you likely still have this heptathlon.csv Assumes PanelLayout functions installed 1. Computing dissimilarity matrices_________________ Dissimilarity matrices and the more stringent special case, distance matrices, have many uses in data transformations and graphics These include Case reduction clustering: an be viewed in this context Dimension reduction multidimensional scaling typically into a low-dimensional Euclidean space principal components non-linear dimension reduction using singular values decomposition methods The notion of dissimilarity is general and can be applied to cases, variables, nodes of a graph, documents, networks, molecular structures, data sets, and so on. It is pretty easy to come up with In most cases the dissimilarity matrix is symmetric with zeros on the diagonal. Cases are not dissimilar from themselves. There are many ways to produce dissimilarity matrices and some follow a minimal set of rules. With scaled data in hand, there still are many ways to compute a dissimilarity matrix from multivariate data. A common approach to comparing cases is to use Euclidean distance. A metric or distance measure for cases i, j and k obeys the following four rules. d(i, j) => 0 : non-negativity d(i, j) = 0 if and only if i = j : identity of indiscernibles. d(i, j) = d(j, i) : symmetry d(i, k) = d(i, j) + d(j, k) : triangle inequality Here it is understood that i, j, and k stand for the multivariate values associated with cases i, j and k. The rules for dissimilarities may drop the triangle inequality, but usually retain the other three. 1.1 Euclidean distance The Euclidean distance between to cases i1 and i2 is the square root of the sum of squares of the difference between variables. As an alternative one my use the pth root of the sum of the pth powers of the differences between variables where p is a even number. If the absolute values of difference are used, p does not need to be an even number. 1.2 Similarities, transformation and correlation In some cases it is convenient to start with an assessment of similarity between pairs of cases. Common similarity rules are s(i,j) = s(j,i) s(i,j) <= s(i,i) Sometimes there is a rule 0 <= s(i,j) <= s(i,i) ==1 1.2.1 Transformations from similarity to dissimilarity A common transformation is d(i,j) = sqrt( s(i,i) - 2*s(i,j) + s(j,j) ) In the case that 0 <= s(i,j) <= s(i,i)=s(j,j) ==1 d(i,j) = sqrt(2*(1-s(i,j)) Dropping the sqrt(2) makes no difference d(i,j) = sqrt((1-s(i,j)) The square root will influence the spacings of of the dissimilarities. As far as I know it is thinkable to skip the square root unless one commmitted precise spacings of the similarities. A second transformation is d(i,j) = 1 - s(i,j)/max(s(i,j)) for i!= j 0 i = j Range [0 1] 1.2.2 Correlation as similarity between variables? Some people have used the correlation, c(i,j), of two variables as an index of their similarity. This a little awkward since correlations can be negative. It is possible to use d(i,j) = 1-c(i,j) Range [0 2] A correlation of 1 maps to 0 dissimilarity which seem reasonable However a correlation of -1 map to the large dissimilarity of 2. I find this hard to interpret. I am more comfortable with d(i,j) = 1 - |c(i,j)| or d(i,j) = 1 - c(i,j)**2 For the normal distribution a zero correlation implies independence. A zero correlation maps into the largest distance of 1. Larger absolute correlations indicate such greater depatures independence. Perfect correlation of 1 or - 1 map into a 0 dissimilarity with seem reasonable. The connection between zero correlation and independence does not carry over to other distribution, but we can still of the distances as reflecting a degree of independence. Other choices include include putting the absolution coorelation |c(i,j)| or square correlation in the first similarity to dissimilarity transformation above. 1.3 1 - cos(angle) between two vectors Range [0 2] cos(angle) is the dot product of two unit length vectors. cos(angle) = (x/||x||) . (y/||y||) Note that if the means have been subtracted from x and y that right hand side is the Pearson correlation. As indicated in 1.2.2 would rather use 1 - |cos(angle)| in this context 1.4 Mutual information Below is an implementation of discrete mutual information function based on partitioning two continuous variables to obtain roughly equal counts in the margins. How many partitions should be used? The default I implemented produces a 3 x 3 table. The partitioning could be much finer. The extreme situation (without ties) has only 1's in the the row and column margins. This can be constructed by putting partitions between adjacent values of the sorted x's and between the adjacent sorted y's). I have not studied the use of mutual information for fine partitions nor mutual information theory in the continuous data context. Using fine partitions might be a bad idea. Mutual information is a measure of similarity. The function has a option to return the dissimilarity value: log2(# of classes) - mutual information ## Run mutualInf = function(x,y,nclass=3,dissim=F){ probs = seq(0,1,length=nclass+1) breaks = approx(ppoints(x,a=1),sort(x),probs)$y xDiscrete = cut(x,breaks,include.lowest=T,labels=F) breaks = approx(ppoints(y,a=1),sort(y),probs)$y yDiscrete = cut(y,breaks,include.lowest=T,labels=F) cnt = table(xDiscrete,yDiscrete) cellProb =cnt/sum(cnt) rowProb = apply(cellProb,1,sum) colProb = apply(cellProb,2,sum) MutualInfMat = cellProb*logb(cellProb/rowProb%o%colProb,base=2) MutualInfMat[cellProb==0] = 0 ans = sum(MutualInfMat) if(dissim) ans = logb(nclass,base=2)-ans return(ans) } x = rnorm(100) y = 2*x + rnorm(100,sd=.2) # adds some contamination # the mutual information goes down # as the contamination increases mutualInf(x,y) mutualInf(x,x) mutualInf(x,y,dissim=T) mutualInf(x,x,dissim=T) # The computed maximum mutual information # does not quite achieve the theoretical value. # I would put 0's in the diagonal elements of # the dissimilarity matrix. ##End 1.5 Manhattan (city block) distance points on a grid 1.6 Jaccard coefficent for binary data The sum of component wise "ands" divided by the sum of component wise "ors" In other words matching zeros are ignored. (1,0,1,0) (1,0,0,1) yields 1/3 not 1/4 In my study of bird species prevalence at many thousands of locations for 600 species, it became clear that the conventional correlation for many species pairs was high because there were so many places where neither of two species were present. It made sense to restricting attention to the locations where at least one of the species was present. The Jaccard coefficent restricts attention to locations were at least 1 bit is non-zero. 1.7 BLended APples and ORanges dissimilarity Blapor? This is a made up name [:-) Cases can be dissimilar in attributes, location and time Ignoring time we could compute dissimilarity as a weighted average of two distances one based on attributes and one based on geospatial location. Note that great arc distance is often more appropriate than Euclidean distance in a map projection This may seem strange and it is uncommon, but dissimilarity values do not have to obey the the triangular inequality. 1.8 Assessing the dissimilarity of different shapes Assessing shape similar (for examples skull shapes) gets interesting. Matching landmarks can be hard. See the Statistical Theory of Shape By Christopher Small 1996, Springer 1.9 More realistic dissimilarities If dissimilarities are viewed as distances between points, they may be made more realistic by considering path constraints between the points. The dissimilarity of concentration measurements in lake may increase with distance through water and opposed to the crow distance over a body of land. If my concern is how long to it takes to drive to work, the number of stop lights can seem as relevant as the distance and the speed limit. 1.10 Nonlinear dimension reduction methods Since 2000, several papers have addressed the ideas of geodesic paths (shortest distance in curved space) and nonlinear dimension reduction. References for 10) Belkin, M. and P. Niyogi. 2003. "Laplacian Eigenmaps for Dimensionality Reduction and Data Representation," Neural Computation, June 2003; 15 (6):1373-1396. Donoho, D. L. and C. Grimes. 2003. "Hessian Eigenmaps: Locally Linear Embedding Techniques for High-dimensional Data, Proceeding of the National Academy of Sciences U S A. 2003 May 13; 100(10): 5591:5596. S. Lafon. 2004. Diffusion maps and geometric harmonics, Ph.D. dissertation, Yale University May 2004. Roweis, S. and L. Saul. 2000. "Nonlinear dimensionality reduction by locally linear embedding," Science, 290(5500), 2323:2326. Tenenbaum, J. B. , V. de Silva, and J. C. Langford. 2000. "A global geometric framework for nonlinear dimensionality reduction," Science 290(5500), 2319:2323. 1.11 Distance between distributions and multivariate distribuitons In statistics there are several way to assess dstance between continous univariate distribution. In earth sciences and many complex application areas we are often be interested in the distance between multivarariate distributions. In class I briefly indicated an expected distance calculation for grid cells on the earth used by Amy Braverman, one of my NASA collaborators. Each grid cell had a set of summary multivariate vectors, each with values for different parameters and altitudes. Each grid cell could have a different number of summary vectors. I used the expected distance to cluster grid cells. This may not be general interest, but I would be glad to go other this in more detail for students with potential applications. 2. The R distance function dist()======================================= The dist() function takes a matrix or data.frame as input. To save storage it outputs the lower triangular part of distance matrix object of class dist. This vector is suitable for input agglomerative clustering using hclust(). ##Run mat = matrix(rnorm(15),nrow=5,ncol=3) matDist = dist(mat) matDist # lower triangle ## End The printing method for a dist object has options to include the upper triangle and diagonal part of the matrixs ## Run print(matDist,upper=T) print(matDist,diag=T,upper=T) ##End If the full matrix is desired for computations ## Run matDistFull = as.matrix(matDist) matDistFull ## End The dist() argument allows specification of several methods for computing distances. This include "euclidean", "maximum", "manhattan", "binary", and Minkowski and a few others. Binary is the same as the Jaccard coefficient as described above, treating all positive values as 1's. If the desired dissimilarity computations are not available, it is usually straight forward to write a function. Obtaining speed and having storage can be problematic for very large data sets. The genes file below may have been used in a previous assignment. If it si not in you workspace can access it from the class schedule. ## Run genes = read.csv(file='gemIdGene.csv',row.names=1) names(genes) genes5 = genes[1:5,1:9] # five genes, times series length 9 geneDist = dist(genes5,method="euc") geneDist round(as.matrix(geneDist),2) ##End 3. Four general approaches to clustering Clustering methods put cases with similar values in the same cluster and cases with dissimilar values in different clusters. In the data mining community clustering often is referred to as unsupervised classification. In supervised classification we have a training set that includes data for cases and their class membership. In unsupervised classification we just have the data and no gold standard from science or experts the class membership for a a training sets. As I sometimes say in class, a new clustering algorithm in created every day. What I mean is that there are lots of clustering algorithms and what works the best can depend on the data set. We are in the position of not knowing what works best. 3.1 Four broad approaches to clustering are 1) agglomorative clustering, 2) recursive partitioning and 3) variants on k-means 4) model based clustering 1)Agglomorative clustering start with individual cases and puts cases and clusters of cases together to make large clusters. Different algorithms have different criteria for merging clusters. Many of the agglomerative cluster algorithms require use of a dissimilarity matrix so are of order O(n^2) where n is the number of cases. There were popular historical but seen less use as data sets have include more and more cases. 2) Recursive partitioning starts with the cases in one large cluster and partitions clusters into smaller clusters. Different algorithms had different criteria for partitioning clusters. 1&2) Dendroograph and color matrices Both agglomerative and recursive partitioning clusters conceptually produce hierarchical structures can often be represented binary trees. The plots of the clustering trees are often called dendrographs. A fairly common plot aligns the leaf nodes of a case dendrogram with the rows (case) of a color matrices encoding the rows and column of data sets. Variables can be clustered as well, so leaf nodes of a variables dendrograms can be a aligned with the column of a color matrix encoding a data matrix. The plot is sometime called a heat map. I many circumstance it is advantageous to show the transpose and the result is still call a heat map. 3) K-means clusters start with means for k-prototype clusters. It the uses a criterion to assign each case to the "nearest" cluster and updates the case means. The protypes many simply be the the mutlivariat values for the first k cases. The prototypes may be the resulting k-means from previously running k-means. The first k cases many not get the algorithm off to a good start, but the hope the ending mean has drifted to a decent starting place. K-means is O(n) so often used for vary large data set. It tends to produce round cluster which is not always desirable. 4) Model based clustering makes assumption the covariance structure of the clusters. Are all the clusterss to be circular or elliptical or some of both? Are all clusters to be roughly the same size or difference sizes in terms of the domain covered. The modeling specification indicate this covariance structure constraints. 3.2 However many clusters should there be? If we are dealing with human cognition a nice number would be four. However we often apply clustering the early stages of science when we don't usually have good basis to make a decision. We know nature is often complex and no respector of human cogntiion. One approach defines a criterion that can be plotted as function of the number of clusters. When high values are bad, the rule of thumb looks for number with big drop followed by diminish returning for using more clusters. Of course this pattern may not occur. Methods for deciding on the number of clusters proliferate. Below are some approach related comments 1) Agglomerative # of clusters After building a tree there is an option to select either the number of clusters desired or a "distance" between clusters beyond which the clusters will not be merged. Examination of the cluster differences may help in making a selection. 2) Recursive partitioning # of clusters One approach conducts hypothesis test to pick the more statistically significant split, adjustes for multliple comparisons and stop partitioning when there no more statistically significant splits. Random forests build on the idea recursive partitioning and can address clustering. This not discuss here but typically in conjugation with a lecture on supervised classification. 3) k-means typically specifies the number of clusters in advance Often k-means will by calculate for several different numbers of clusters to provide a based for comparison use a cluster criterion as indicated above. One variant of k-means specifies a large number of clusters such as 2^10 = 1024 and keeps tract of the means and covariances. Then procedure pools close clusters to reduce the number of clusters. These has been used on some very big data sets. Another variant is called entropy constrained vector quantization ECQV. Here the analysts specifics a maximum tolerated number of clusters (sometimes after look as clustering results for parts of the data). If fewer clusters statisfy and entropy panel then few are return. The class example on global multivariate multialtude atmospheric data compression used ECQV. The number of cluster summarizes for each 5 x 5 degree grid cell on the varied. This was applied to terabytes of data. The first 18 principal components of 35 variables captured almost all the variability so were used in clustering. The given the cluster membership for satellite footprint the last processing step computd the means using the geophysical parameters. My impression is that some variant of k-means is almost always used for very large data sets because most of the alternative methods require too much in the way of computer resource (time, storage, or both). 4) I don't have special number of cluster comments for model based clustering. In general I like the idea of incorporationg additional knowledge into the process whenever such is available. 3.3 Cluster evaluation Simulation studies can construct clusters and compare error rates for different clustering methods. In the absence of "truth," clustering methods can still be compared and assess in various way. In collaborate research at GMU I provided several ways of looking gene expression time series in the hope of encouraging graphical evaluation or criticism. Some appear a newsletter article "Templates for looking at Gene Expression Clustering" by Carr, Michaels, Somogyi, and Micheals [1997]. that is available in V81.pdf from the class web site. After scaling the time series and clustering the result using Euclidean distance, I displayed the time series residuals from the mean times series of each cluster in a two way layout of panels. There were different rows for different gene families and different columns for the different clusters. We had chosen to use 6 clusters so there were six columns. The objective was to compare the clusters with exisiting scientific knowledge. A natural question is: "is there some kind of match between gene family membership and cluster membership?" This translates into a visual pattern finding query of looking for panels with small residual series. A related plot sometime give in a class assignment compare two choices for dissimilarity measures using this data sets. We had used both Euclidean distance and mutual information as a basis for clustering. The comparative plots was row and column panel layout that overplotted the time series for the genes. For examples time series for genes in the 1st Euclidean distance cluster and 2nd mutual information cluster appeared in the (1st row, 2nd column) panel and so on. A natural question is, "Are some of the clusters using the different dissimilarity measure basically the same?" This translates into visual query scan for a panel with limited variation in the a panel and relatively few genes in the remaining panels of the same row and the same column. Another kind of scan looks for =big variation in any panel that suggests the clusters were not very good (tight). A second interesting plot in the paper was the side by side 3D scatterplot. This showed the first 3 principal components for the time series. It distinguised the 6 Euclidean distance clusters with 6 colors. The plots also connected the points in each cluster with minimal spanning tree lines. Some clusters seem tight and some looked sprinkled across a large part of the 3D domains. I chose to minimize overplotting so associated the 1st, 2nd and 3rd principal components with x(width), y(height)and z (depth) in the side by side stereo plot. The axes lengths conveyed the ranges of the three principal components. An interesting choice to consider is the choice spanning tree used to to connect the points. If the spanning used is constructed using the first 3 principal components the plot appears consistent. Connecting points based on the spanning tree in the full dimensions of the data can bring bring out possible coordinate anomalies. In Crystalvision I would pick 3 time series coordinates (early, middle and late rat ages), color the points by cluster membership and draw the full dimension spanning tree lines. Some points and their lines look inconsistent. There are different possibile explanations. For example a value for one of the coordinates used in the plot may be bad but not bad enough to pull the gene out of the cluster. An interesting construction challenge was to draw the color lines so the the correct color was on top in the absence of a z-buffer. I divided the lines into tiny segments and plotted the segments back to front. This removed most of the errors. (A disgruntled student brought down a lab in STII and erase my Splus script file.) A third interesting plot in the paper made it into an GMU science and art exhibit. This suggest a kind of meta analysis that I call coherence analysis. Perhaps this kind of higher order analysis has been explored but as best I know only the coauthors have also thought about the broader implications. It is not just a pretty plot. 3.4 Finding subspaces with clusters Some variables may be irrelevant to clusters and different sets of variables may be relevant to different clusters. I encounter this topic when collaborating with Alan MacEachren (Penn State) and interacting with some of his students. Later I learned Carlotta Domeniconi and Daniel Barbara in IT&E at GMU had been developing algorithm to consider different local subspaces in constructing clusters. If someone is really interest I suggest contacint Dr. Domeniconi. 4. A look at agglomerative clustering with hclust() Section 4.1 provide an introduction to clustering based on old hclust() algorithm for clustering. I still use this function on occasion and the data structures can be interesting. 4.1 The hclust() options and output_________________________ The hclust() clustering methods are "ward" , "single","complete","average", "mcquitty","median",or "centroid" The methods are the same for joining two singleton items. They select the pair of items with the smallest dissimilarity. The methods also join the two groups with the smallest dissimilarity but differ in how this dissimilarity is computed. Note the one of the "groups" can have just one case.. Single: For the single link method, the dissimilarity between groups is is the smallest dissimilarity considering pairs of items, one from each group. This approach makes it easier to prove some theoretical properties and will be illustrated below because it is easy to follow. However, this approach tends to produce long chains that join singletons to groups one at a time. Applied researcher have not found the results appealing. Compact: Little used Average: The group dissimilarity for the average method is the average of dissimilarities between pairs using one item for each of the two groups. Applied researchers generally prefer this over single link and compact algorithms. ##Run clustAns = hclust(geneDist,method='single') ##End The output is a list with three elements, a merge matrix, a merge height vector, and the item plotting order for cluster tree graphics historically called dendrograms. ##Run clustAns$merge ##End [,1] [,2] [1,] -3 -4 [2,] -1 -2 [3,] 1 2 [4,] -5 3 !!! note there are 5 items and 4 joins Negative values in the two column matrix indicate the original points Positive values in the matrix indicate the row of the matrix that stands for a cluster Cluster 1 has points 3 and 4 Cluster 2 has points 1 and 2 Cluster 3 is composed of clusters 1 and 2 Cluster 4 is composed of point 5 and cluster 3. In general n items and n-1 joins ##Run round(clustAns$height,2) ##End 0.48 0.49 0.52 0.71 The above vector contains the four dissimilarities between merged points or clusters being joined ## Run clustAns$order ## End [1] 5 3 4 1 2 The above gives the ordering of cases to produce a cluster tree or dendrogram ## Run plclust(clustAns) ## End Why not transpose this dendrgram so the labels are easier to read. 4.2 Stepping through single link clustering____________________ ## Run dissim = as.matrix(geneDist) round(dissim,2) ## End [,1] [,2] [,3] [,4] [,5] [1,] 0.00 0.49 0.52 0.73 0.94 [2,] 0.49 0.00 0.68 0.55 0.97 [3,] 0.52 0.68 0.00 0.48 0.72 [4,] 0.73 0.55 0.48 0.00 0.71 [5,] 0.94 0.97 0.72 0.71 0.00 To start, each item is considered a group. The closest two "groups" are items 3 and 4 with a dissimilarity of .48 The first line of ans.con$merge indicates we join these in a group G1 The first element in ans.con$height is the merge height, .48. Now we update the dissimilarity matrix. The single link group dissimilarity is the smallest dissimilarity between the pairs of elements, one from each group. G1 = c(3,4) Dis 1 to G1 = min( dis(1,3),dis(1,4) ) =.52 Dis 2 to G1 = min( dis(2,3),dis(2,4) ) =.55 Dis 5 to G1 = min( dis(5,3),dis(5,4) ) =.71 The updated dissimilarity matrix is [,1] [,2] [,G1] [,5] [1,] 0.00 0.49 0.52 0.94 [2,] 0.49 0.00 0.55 0.97 [G1,] 0.52 0.55 0.00 0.71 [5,] 0.94 0.97 0.71 0.00 By inspecting the matrix or by looking at the second line of ans.con$merge we see that next closest groups are singletons 1 and 2 with a dissimilarity of .49. This is the second element in ans.con$height We join these in a group G2 and update the dissimilarity matrix. G2 = c(1,2) Dis G1 to G2 = min( dis(G1,1),dis(G1,2) ) =.52 Dis 5 to G2 = min( dis(5,1), dis(5,2) ) =.94 The new dissimilarity matrix is [,G2] [,G1] [,5] [G2,] 0.00 0.52 0.94 #.52 is the smallest values [G1,] 0.52 0.00 0.71 [5,] 0.94 0.71 0.00 The closest groups are now G1 and G2 with a dissimilarity of .52 We join these created group G3 and update the dissimilarity matrix G3 = c(G1,G2) Dis (5,G3) = min(dis(5,G1),dis(5,G2))=.71 [,G3] [,5] [G3,] 0.00 0.71 [5,] 0.71 0.00 The only groups are now G3 and 5 with a dissimilarity of .71 We join these in group G4 # Check round(ans.con$height,2) [1] 0.48 0.49 0.52 0.71 4.3 Merge order for dendrograms and linear orderings of multivariate cases_________________ #clustAns$merge # # [,1] [,2] # [1,] -3 -4 # [2,] -1 -2 # [3,] 1 2 # [4,] -5 3 Note three construction details. 1) If two single item are merged, the smaller absolute row subscript of the data matrix appears first. 2) If two groups are merged the smaller group number appears first. Since the smaller group number cluster was merged first it was "tighter" cluster. 3) If a single item and a group are merged the single item appears first. # Following the recipe in ans.con$merge yields ##Run G1 = c(3,4) G2 = c(1,2) G3 = c(G1,G2) G4 = c(5,G3) ##End G4 [1] 5 3 4 1 2 This is the same as ans.con$order #[1] 5 3 4 1 2 As indicated above this is a valid case ordering for producing a dendrogram. Valid means that connecting merged groups at the merge heights does not involve crossing lines when using vertical and horizontal lines. See the example below ## Run windows() plclust(clustAns,hang=-1) ## End Note the arbitrary choices made in the construction details. Various flips of the pairs would still produce valid dendrograms. The ordering 5 4 3 1 2 would work as well as 5 3 4 1 2 as well as 5 3 4 2 1 as well as 3 4 2 1 5 The ordering of the genes in the original data matrix was arbitrary but nonetheless influences the default dendrogram ordering. In some situtaions switching the pairings above reduces the sum of consecutive dissimilarities in the linear ordering of a dendrogram. # A two column matrix of subscript pairs can be used to # extract matrix elements. ##Run # original order subs1Mat = cbind(G4[-length(G4)],G4[-1]) subs1Mat # #integer matrix: 4 rows, 2 columns. # [,1] [,2] #[1,] 5 3 #[2,] 3 4 #[3,] 4 1 #[4,] 1 2 # new order subs2Mat = matrix(c(5,4, 4,3, 3,1, 1,2),ncol=2,byrow=T) # Remember dissim is a symmetric matrix from above dissim[subs1Mat] # [1] .7193747 0.4795832 0.7312318 0.4904080 sum(dissim[subs1Mat]) # Sum of dissimilarities for adjacent genes # 2.42 sum(dissim[subs2Mat]) # Sum of dissimilarities for adjacent genes # 2.20 ##End Here switching the order of the 3rd and 4th genes produces smaller adjacent dissimilarities in the linear sequence. A dendrogram ordering provides a linear ordering of multivariate cases. The task of linearly ordering multivariate observations (see also seriation) is recurrent task in quantitative graphics. My experience hasled me to prefer the minimal spanning tree breadth traversal ordering for this purpose. I have been using this to order variables in a correlation matrix from the mid 1980s not that anyone noticed. See the Carr and Olsen 1996 newsletter article "Simplifying Visual Appearance By Sorting: An Example using 159 AVHRR Classes. In more recent times Stan Young suggested that I try using the first left eigenvector of a singular value decomposition to rearrange the rows and first right eigenvector to order columns. In some instances prefer this but more often prefer the minimal spanning tree traversal order. There is an Splus function clorder() that can be used to obtain a alternative dendrogram orderings. This could be used in an attempt to get closer to the minimal spanning tree traversal. There is a spanning tree package for R. As far as I know the breadth minimal spanning tree traversal algorithm described Friedman, J. H. and Rafsky, L. C. (1981). "Graphics for the Multivariate Two-Sample Problem." Journal of the American Statistical Association 76, 277-287. is still not available for R. I have started developing such an algorithm but have not yet finished it. 5. Dendrograms and cluster tree variations__________________________ ##Run nams = row.names(genes5) plclust(clustAns,hang=-1,labels=nams) nams ##End [1] "actin" "SOD" "CCO1" "CCO2" "SC1" Note that is not so easy to read the vertical labels. Also observe that the 5th item appears first in the tree based on ans.con$order In looking at the plots, the distance from the leaves to the first merge can seem to consume much space. An alternative is to draw cases at fixed distance down from the merge height. This is the hang parameter. By default it is .1 of the plot height. ##Run plclust(clustAns,labels=nams) ##End If the merger order is more important than the merge height, the integer order can be on the y axis. ##Run plclust(clustAns,labels=nams,unit=T) ##End ## Skip since Splus options were not yet working in R the last time I checked__________________ To collapse this for large trees there is a leveling option plclust(clustAns,labels=nams,unit=T,level=T) For larger trees the merge details below a certain merge height may not be of interest. The parameter hmin controls the hiding of merge details. Here are two variations plclust(clustAns,labels=nams,hmin=.6) par(mai=c(1,1.5,1,.5)) plclust(ans.con,labels=nams,hmin=.6,xlab='Genes',ylab='Merge height',main= 'Merge Tree') # End Skip Splus options________________________________________ Horizontal labels are nice. This motivates a horizontal tree orientation. Also smooth lines such as parabolic lines are easier to follow than following lines around corners. Unfortunately they can cross on occation. Get my hypertree function form the web site. (When wrote this the function I had planned to use hyperbolic curves hence the label.) ##Run source('hypertree.r') windows() hypertree(clustAns,lab=nams,xsh=.005,padx=0.02) ##End Argument include xsh: controls how close the text is to the tree padx: controls the space for the labels drop: controls the drop height (same purpose as the hang parameter above) 6. A dendrogram for a subset of gene expression data================== Each row gives multivariate values for a gene. The genes are a tiny set of those controlling the development of the rat spinal chord. The first nine columns are a time series. These are scaled values describing the level of gene expression (mRNA production) at different times during the life of a rat. The time label interpretation: E11 means gestation day 11. P0 means day of birth p14 means 14 days after birth A means adult. 6.1 Some general background on genes and gene expression The expression levels reflect the amount of mRNA produced by the genes. mRNA production (transcription) is a part of the protein production process. The role of some proteins is to regulate gene expression levels. The proteins dock upstream of a gene's coding region and up or down regulate the coding that produces mRNA. Several kinds of proteins call transcription factors can be involved in the transcription regulation of a gene. Sometimes regulation applies to a group of genes collocated on a chromosome. When mutation moves some genes to another location on a chromosome or to another chromosome, the regulatory process performs poorly or fails. Amazingly, not all of such major mutations are fatal. There are living human beings with groups of genes on the wrong chromosome! Many of these people have serious problems. While every cell has the same DNA, expression levels and regulation differ for different tissue types in the body. A cell of the eye responds differently than a cell in the toe or a cancerous cell in liver. Genes are part of metabolic pathways that are regulated by both external events and other genes. Understanding the communication process in the context of multiple biological systems (such as metabolic pathways) remains a huge scientific challenge. Looking the the broad scope of communication is call the Systems Approach 6.2 More on the rat data In general each value in the data set was based on expression levels from three rats. The individual values and thea summary variability measure are not available here. The experiment sacrificed the rats to obtain data. Values obtained at different times came from different rats. The biological variability from animal to animal can be huge. 6.3 Scaling the expression values_____________________________ If all the multivariate observations are measured using the same units, it may not be necessary to scale the data. Often a first step is to "scale" the data into unitless numbers. For example dividing by values of each variable by the variable standard deviation will work. The nine values for each gene have been transformed so the maximum value is 1. Some genes may start near the maximum production level and stay at the level for the life of the animal. This motivates scaling by just dividing by the maximum value. This transformation is much less common than the transformation that scales variables into the interval [0 1] by first subtracting the minimum value. 6.4 Looking at a bit of data # The command below obtains the gene expression data for 5 genes at 9 times. ##Run dat = genes[1:5,1:9] # redirect the text to a file in the workspace sink('dat.txt') dat sink() # point the output back to your monitor ndat =as.matrix(dat) # convert the data frame to a matrix # for use in dist() below nams = row.names(dat) # store the gene names for each row ##End 6.5 Creating new variables_________________________ The investigators thought that changes between times between the ages should be emphasized. They could have expressed the changes of slopes that took into the time difference. However they chose to just augment the data with differences. ##Run augMat = t(apply(genes[,1:9],1,diff)) # conversion of data.frame to matrix geneMat = cbind(genes[,1:9],augMat) ##End 6.6 A dendrogram for a larger subset of gene data ##Run gene30Mat = geneMat[1:30,] gene30Nams = row.names(gene30Mat) gene30Dis = dist(as.matrix(gene30Mat)) gene30Clus = hclust(gene30Dis,method='ward') windows() hypertree(gene30Clus,lab=gene30Nams,xsh=.005,padx=0.02,cex=.55) ##End # Note if the "ave" clustering is used above it reveals that # the hypertree algorithm can have crossing lines 6.7 Comment on graphics goals Clearly it is going to be difficult to study cluster trees from leaf nodes to the top for very large data sets. There are n leaves and n-1 joins. Sometimes it is useful to look clusters at the top of the tree to see their relationships. In class I show a variant of the hypertree graphics limited to the top of the tree. The example shows the tree top for the Stanford yeast gene dataset with over 6000 genes. The example includes that overplotted expression times series (or profiles) for the 30 cluster "leaf" nodes. One cluster, the non-active cluster, involved over 5000 genes but expression pattern was clear. I can make the additional scripts available to interested students. Goals include cluster description, understanding and criticism. There are so many things that have been done and additional possibilities. Suppose there some bad measurements. If there is some redundancy in the variables, we could compare clusters against those based on leaving out each variable in turn. This could point to cases with an inconsistent observations. One interesting approach to studying clusters uses a space filling layout to plot cases as points in a scatterplot so that rectangles can enclose clusters at any level of clustering. (I gave the space filling layout construction as an SCS computational exam question, and this was later used in the students dissertation on document clustering). This plot combined with a scatterplot matrix of the variables provides a frame work for linked brushing. More recently I develop a hexagon cluster graph to use for brushing and displaying glyphs. This was illustrate in class in the context of visualization clusters 5 x 5 degree cells on the earth talk about It used to be that statisticians would expect the number of cases to exceed the number of variables by at least a factor of 3. Today there may be many more variables than cases. For example document clustering may be based on word counts for unique words numbering in the several tens of thousands. In some cases counts for each unique word represents a variable. In other cases counts for each unique pair of words (bigrams) represent a variable. There will be a more examples related to clustering in class and possibly in latter exercises. However there are additional graphics topic that need attention Right now we move on to defining clusters 7. Defining clusters using cutree() The cluster tree shows all items merged into one cluster. We create clusters by not making the last merges. There are two standard ways to prevent the last merges. 7.1 Specify the number of clusters. Specifying k clusters prevents the last k-1 merge. ##Run gene30Grps4 = cutree(gene30Clus,k=4) gene30Grps4 # # actin SOD CCO1 CCO2 SC1 SC2 SC6 SC7 # 1 1 1 1 2 1 3 1 # # DD63.2 cyclin_A cyclin_B H2AZ statin cjun cfos Brm # 1 1 2 1 1 4 4 2 # TCP CRAF IP3R1 IP3R2 IP3R3 Ins1 Ins2 IGF_I # 1 1 1 1 2 3 1 1 # IGF_II InsR IGFR1 IGFR2 NGF NT3 # 2 1 2 2 4 2 table(gene30Grps4) # There are 17 genes in the first group #gene30Grps4 # 1 2 3 4 #17 8 2 3 gene30Grps4[gene30Clus$order] # dendrogram order ##End 7.2 Define clusters by merge height The second method specifies the maximum allowable merge height. Any group that would have merged at greater heights (dissimilarities) is retained as a group. ##Run gene30ClusterId = cutree(gene30Clus,h=1.2) # merge up to a height of 1.2 gene30ClusterId # actin SOD CCO1 CCO2 SC1 SC2 SC6 SC7 # 1 1 2 2 3 4 5 4 # DD63.2 cyclin_A cyclin_B H2AZ statin cjun cfos Brm # 1 6 7 6 8 9 9 3 # TCP CRAF IP3R1 IP3R2 IP3R3 Ins1 Ins2 IGF_I # 6 10 1 1 11 12 2 6 # IGF_II InsR IGFR1 IGFR2 NGF NT3 # 3 10 7 7 9 13 ## End At this height there are still 13 clusters ##Run table(gene30ClusterId) ## End gene30ClusterId 1 2 3 4 5 6 7 8 9 10 11 12 13 5 3 3 2 1 4 3 1 3 2 1 1 1 #7.3 subtree() #Not available in R # #The subtree() function returns the smallest subtree that #contains a give set cases. It can be applied to focus on a #cluster by indicated the correct case subscripts. # ###Run #clus1 = seq(along=gene30ClusterId)[gene30ClusterId==1] #clus1Nam = gene30Nams[clus1] #gene30Cluster1Clus = subtree(gene30Clus,clus1) # #gene30Cluster1Clus$merge ###End # # [,1] [,2] #[1,] -16 -25 #[2,] -5 1 #[3,] -27 -28 #[4,] -11 2 #[5,] 3 4 #[6,] -21 5 # #This will work with plclust but not with the hypertree #function that expects consecutive case numbers starting #from 1. 8. Heat maps and Color Matrices, Cluster Visualization Link I think of color matrixes and heat maps as the same thing. Colors matrices can show the data values, scaled data values, correlations, dissimilarities, etc. Heatmaps sometimes show data values with dendrograms on two sides of the matrix. Perhaps that is the distinction. Another option is to use color along two sides to show dissimilarities for adjacent cases and adjacent variables 8.1 A R example heatmap ## Run require(graphics) x = as.matrix(mtcars) rc = rainbow(nrow(x), start=0, end=.3) cc = rainbow(ncol(x), start=0, end=.3) hv = heatmap(x, col = cm.colors(256), scale="column", RowSideColors = rc, ColSideColors = cc, margin=c(5,10), xlab = "specification variables", ylab= "Car Models", main = "heatmap(, ..., scale = \"column\")") str(hv) # the two re-ordering index vectors ##End This could used better colors and be refined. See the R help on heatmap for more examples and explanation. 8.2 Optional: interactive heat map visualization Those with windows might want to give the software with access below a try. Go to http://www.cs.umd.edu/hcil/hce/ Register and download the software. I suggest the HCE3.5 version. You might also want down load the manual. When you run HCE you might start by selecting the cereal example. I usually transform the columns to z scores. You can cluster both rows and columns. I will show some of the options in class. NCI was thinking about and perhaps used the software to identify sister counties to used on comparative cancer studies. 9. Clustering other data sets Sometimes it is suggestive to cluster data when the cases in a data set has already been classified. Clustering can be useful in various ways but wild enthusiasm for clustering may lead to some to strange applications. This assignment closes by clustering point in a hypersurface, a strange application. 9.1 Fisher's iris data The data, called iris, already in R. It is stored as a data.frame giving 4 measurements and species type for 150 flowers. The measurements in centimeters are of sepal length and width, and of petal length and width. There are 50 flowers for each of 3 species of iris. The species are Setosa, Versicolor, and Virginica. The species is stored as a factor R. A. Fisher, "The Use of Multiple Measurements in Taxonomic Problems", Annals of Eugenics, 7, Part II, 1936, pp. 179-188. ##Run iris irisDis = dist(iris[,1:4],method="e") irisClus = hclust(irisDis,method="ward") irisClass = cutree(irisClus,k=3) table(iris[,5],irisClass) ## End irisClass 1 2 3 setosa 50 0 0 versicolor 0 50 0 virginica 0 14 36 The cluster Class numbers 1-3 are arbitrary and could be switched to create more of a diagonal matrix. The about matrix is fine as is. The 1st cluster include exactly setosa flower The 2nd cluster picks up all versicolor flowers The 3rd cluster up a little over half of the virginica In other words the clustering has matched the scientific knowledge quite well for two of the three clusters Other clustering methods may do better or worse. ## Run windows() plclust(irisClus) ##End Note the there can bad overplotting in the plclust algorithm. The HCE software shown in class has many merits for example handling larger data sets. 9.2 Clusters and statistics for athletes in the heptathon data Get the data from the class web site ##Run hepDat = heptathlon[,-c(1,ncol(heptathlon))] # remove names and final score heptathlonDis = dist(hepDat,method='e') heptathlonClust =hclust(heptathlonDis) windows() plclust(heptathlonClust,labels=heptathlon[,1]) # A look at means and covariances for 4 clusters heptathlonClustId = cutree(heptathlonClust,k=4) table(heptathlonClustId) grpLab = c('G1','G2','G3','G4') varLab = colnames(hepDat) nCol = ncol(hepDat) clusMeans = matrix(0,nrow=4,ncol=nCol) dimnames(clusMeans) = list(grpLab,varLab) clusCovars = array(data=0,dim=c(nCol,nCol,4)) dimnames(clusCovars) = list(varLab,varLab,grpLab) for(i in 1:4){ mat = hepDat[heptathlonClustId==i,] clusMeans[i,]= apply(mat,2,mean) clusCovars[,,i] = var(mat) } # cluster means round(clusMeans,2) # cluster covariances in an array clusCovars[,,1] clusCovars[,,2] clusCovars[,,3] clusCovars[,,4] ## End !!! Note that one cluster has only one case so there is no covariance matrix. 9.3 Hypersurface generation, clustering, and graphics 9.3.1 Hypersurface construction A surface can be a described as three functions of 2 parameters: u1 = f1(x,y) u2 = f2(x,y) u3 = f3(x,y) A hypersurface can be described as four functions of 3 parameters u1 = f1(x,y,z) u2 = f2(x,y,z) u3 = f3(x,y,z) u4 = f4(x,y,z) To construct random points on a hypersurface we 1) concoct four continuous functions. 2) Repeat n times Generate random values for the three parameters Use the functions to calculate four values given the three parameters Below ##Run This is organized to avoid looping # Generate 3 parameters x = runif(200,-.2,.8) y = runif(200,-.2,.8) z = runif(200,-.2,.8) # Generate 4 values using made up quadractic functions u1 = 3*x^2-2.5*x*y-.2*z u2 = -.3*x^2 + 2*y^2 +.5*z^2 + 4*x*y -1.2*x*z -.2* y*z + + x -1.1*y +1.3*z +1 u3 = y-3*z^2 u4 = 7*y^2 -3*x*y + 4*x*z -.09*z^2 hyperMat = cbind(u1,u2,u3,u4) ##End 9.3.2 Clustering a hypersurface and a scatterplot matrix ## Run hyperMatDist = dist(hyperMat) hyperMatClus = hclust(hyperMatDist,method="average") hyperMatClusId = cutree(hyperMatClus,k=4) pairs(hyperMat,col = hyperMatClusId,pch=16,gap=0) ## End 9.3.3 Define function for rotating ray glyphs We should discuss the scaling in the function in class. Note the function support appending variable means for cluster to the bottom of the data matrix. Argments groups: Give the cluster id for each case it also give the cluster id for group means if they are added to the data there are cluster means all extra: How many rows were added to the matrix. Typically this is the number of clusters pch: Can have two values, one for data and one for added points cex: Can have two values, one for data and one for added points A different color index is used for each group. The colors depend on the default palette ## Run rayRotate=function(mat,nviews=500,groups=1,extra=0, pch=20, dotSize=.45,rayLength=.03,wait=.03){ # rayLength is the fraction of the half range # Scale in [-1 1], subtract the midrange as a start matRange = apply(mat,2,range) matScaledOne = scale(mat,center=(matRange[1,]+matRange[2,])/2, scale=diff(matRange)/2) apply(matScaledOne,2,range) # check x = matScaledOne[,1] y = matScaledOne[,2] z = matScaledOne[,3] w = matScaledOne[,4] d = sqrt((matScaledOne[,1:3]^2) %*% c(1,1,1)) # largest value from origin of 3D rotation d = max(d) plim = c(-d,d) dx = rayLength*d*cos(w*pi/2) dy = rayLength*d*sin(w*pi/2) windows(width=8,height=8.3) par(pty='s',mai=c(.02,.02,.02,.02)) nviews = 500 n=nrow(mat) good = c(rep(T,n-extra),rep(F,extra)) left = !good for (i in 1:nviews){ angle = (.02*i) %% (2*pi) xnew = x*cos(angle)+z*sin(angle) mark = proc.time()[3] plot(plim,plim,type='n',axes=F,xlab='',ylab='') points(xnew[good],y[good],pch=pch[1],col=groups[good],cex=dotSize[1]) if(any(left))points(xnew[left],y[left],pch=pch[2],col=groups[left],cex=dotSize[2]) segments(xnew,y,xnew+dx,y+dy,col=groups,lwd=2) waited = proc.time()[3]-mark while(waited < wait)waited = proc.time()[3]-mark } } ##End ## Run clusMeans = matrix(0,nrow=4,ncol=4) for (i in 1:4){ mat = hyperMat[hyperMatClusId==i,] clusMeans[i,] = apply(mat,2,mean) } hyperMatPlus = rbind(hyperMat,clusMeans) groups = c(hyperMatClusId,1:4) extra=4 pch=c(20,20) dotSize=c(1.5,2.5) rayRotate(hyperMatPlus,nview=500,groups=groups,extra=4,pch=pch, dotSize=dotSize) ## End You may conclude from this example that clustering is not very helpful for finding key features of data in some context 10. Kmeans K-means is easy to use although there are many options. Centers= is either the number of clusters (# of cases to be randomly selected as beginning centers from the data) or a matrix of distinct starting centers. iter.max = the number of times to pass through data nstart = the number of times pick with random cases as starting cluster centers. In the example below the algorithm will start at 3 different random sets of starting points. How many clusters to use is a recurrent question. ##Run nclust = 4 kmeanAns = kmeans(hepDat,centers=nclust,iter.max=5,nstart=3) nclust= nrow(kmeanAns$centers) pairs(heptathlon, col = kmeanAns$cluster,pch=16,gap=0) ##End